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Abstract 

We present a dedicated analysis of the influence of excited states on the 
calculation of nucleon matrix elements. This calculation is performed at a 
fixed value of the lattice spacing, volume and pion mass that are typical of 
contemporary lattice computations. We focus on the nucleon axial charge, 
gA, for which we use about 7,500 measurements, and on the average mo- 
mentum of the unpolarized isovector parton distribution, {x) u _ d , for which 
we use about 23,000 measurements. All computations are done employing 
Nf = 2 + 1 + 1 maximally-twisted-mass Wilson fermions and using non- 
perturbatively calculated renormalization factors. Excited state effects are 
shown to be negligible for gA, whereas they lead to an O(10%) downward 
shift for (x) u _ d . 

Keywords: Nucleon Matrix Element, Parton Distribution Function, 
Lattice QCD 



1. Introduction 

Understanding nucleon structure is one of the fundamental goals of lat- 
tice QCD. Such an endeavor is becoming more realistic as present day cal- 
culations are being performed closer to the limit of physical quark masses, 
small lattice spacings and large volumes pQ . Thus a direct comparison of re- 
sults from lattice calculations and experimental measurements is becoming 
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feasible, allowing us to probe QCD as the underlying theory of the strong 
interactions. 

In order to establish that lattice QCD can indeed provide results that 
address this challenge, we are focusing, in this work, on a number of bench- 
mark observables for which QCD is expected to produce the right results. 
However, even with the significant advances of the past few years, there is 
presently an unexpected discrepancy between lattice calculations and ex- 
perimental measurements of physical observables such as the nucleon axial 
charge, gA, the average momentum of the unpolarized isovector parton dis- 
tribution, {x) u _ d , and the charge radius of the nucleon [2H5]. Clearly, these 
discrepancies need an explanation. Naturally, this demands a careful study 
of the systematic effects that play an important role in lattice calculations, 
i.e. lattice artifacts, finite volume effects and non-perturbative renormaliza- 
tion. Indeed, these systematic effects are currently being addressed by many 
lattice QCD collaborations. In this work, we focus on the investigation for 
two key observables, namely gA and {x) u _ d . These quantities are deter- 
mined on the lattice by computing suitable ratios of 3-point and 2-point 
correlation functions that reach a constant plateau value for asymptotically 
large Euclidean time separations. 

There are several studies of gA and (x) u _ d at different values of the lat- 
tice spacing and various volumes [BHS]. For the values of the pion masses 
considered, these results reveal that taking into account lattice spacing and 
finite volume effects is probably not sufficient to reconcile the lattice calcu- 
lations with the current experimental value of gA and the phenomenological 
extractions of (x) u _ d . In addition, there is a recent study indicating that 
the discrepancy persists even for pion masses quite close to the physical 
value [9]- One might still argue that calculations performed at precisely the 
physical point might eliminate these discrepancies but that would require a 
strong pion mass dependence making such an explanation increasingly less 
plausible. 

An issue that is currently under scrutiny within the lattice community 
concerns excited state contamination. These are sub-leading contributions 
in Euclidean time correlation functions that can cause a systematic effect in 
determining the desired nucleon- nucleon matrix element. In order to under- 
stand whether the physical plateau appears at larger Euclidean times than 
used in present day calculations, we need to examine the relevant 3-point 
functions at Euclidean times that are as large as possible. The difficulty as- 
sociated with probing such sub-leading contributions to these ratios is that 
the signal-to-noise ratio decreases exponentially fast with Euclidean time. 
To tackle this problem, we have performed a dedicated high precision cal- 
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culation of the 3-point functions for gA and {x) u _ d . In particular, we have 
used about 7500 measurements for gA and about 23000 for {x) u _ d . With 
such large numbers of measurements, we are able to calculate the correlation 
functions for gA and {x) u _ d for large Euclidean times to sufficient accuracy. 
This allows us to study carefully possible excited state contributions, which 
is the goal of this work. 

This paper is organized as follows. In section [2j we describe briefly how 
the matrix elements for the observables of interest in this work are calculated 
in Euclidean field theory and how the excited state contributions enter the 
calculation. In section [3j we give the details of the calculation and explain 
the open sink method used in this dedicated analysis. In section [4] we present 
our results and in section [5] we summarize and conclude. 

2. Nucleon Matrix Elements and Excited State Contributions 

In order to make the paper self-contained, we introduce the basic defini- 
tions of the quantities studied here. In the following discussion, we assume 
an infinite volume, while in all our practical calculations, periodic or anti- 
periodic boundary conditions are taken as needed. The effects of a finite 
space-time lattice are felt in the standard finite-size effects of the matrix 
elements and the so-called thermal contributions that distort the Euclidean 
time dependence of correlation functions. These systematic effects are both 
addressed by the finite-size studies discussed previously and are henceforth 
ignored. 

The zero momentum nucleon 2-point correlation function on the lattice 
is defined as 



where Jat is a nucleon interpolating field that creates a state on the lattice 
with the same quantum numbers as the nucleon, T is a matrix acting in Dirac 
space, and the sum over the Dirac indices a is implicitly understood. The 
Euclidean time t denotes the separation between the creation time t S ource 
and annihilation time i s i n k of the nucleon and is often referred to as the 
source-sink separation, t = t s ; n k — t source . We use translational invariance to 
set tsource = resulting in t = t s ink- The nucleon interpolating field is 
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where C = ^7072 is the charge conjugation operator. The transfer matrix 
formalism allows us to relate lattice correlation functions to matrix elements 
of operators. Application of the standard methods gives the spectral repre- 
sentation of the 2-point function in terms of the eigenstates of the transfer 
matrix or equivalently the Hamiltonian %. The resulting expression is 

C 2 (t) = ^J«j}f ) e-^ t , 
k 

which in the limit t —> 00, will be dominated by the nucleon ground state, 

limC 2 (i)4 J^J^e- mNt . (1) 

In these expressions, labels the masses in the nucleon channel and 

denotes specifically the nucleon mass. Additionally, we have introduced the 

ik) ~(k) t> 

symbols J^- and for the overlap of the interpolating fields with the k 

eigenstate of %. Strictly speaking, the limit t — > 00 cannot be realized on a 

finite lattice; in practice however, it suffices to take t large enough so that 

the correction coming from the lowest lying excited state can be neglected. 

For the evaluation of nucleon matrix elements that we are interested in, 

we additionally need to calculate 3-point correlation functions. They are 

defined as 

c 3 (t, t') = r' aa , (jN,a>(x, t) 0(y, t') J N , a (6, 0)) , (2) 

where is a local field corresponding to the operator O of interest and 
r' is an appropriately defined matrix acting in Dirac space. We denote by 
t' the insertion time of the operator under consideration. Like the 2-point 
function, there is a spectral representation that can be derived from the 
transfer matrix formalism and reads 

^3 (t, f ) = £ 4 ] Jjj*>e-^<*-*V™**' (j| 6 \k) . 

The asymptotic limit of large Euclidean time again isolates the nucleon 
contribution as follows 

Hm C 3 (t, f ) -»• jWjWe-™* (0| 6 |0) . 
(t—t'),t'-* 00 

The desired nucleon matrix element (0| O |0) is then obtained from the 
asymptotic Euclidean time limit of the ratio of the 3-point and the 2-point 



4 



function 



<O|0|O)= Hm ^MO. (3) 

It is the main goal of this paper to investigate how large t — t' and if should 
be so that the contribution of the lowest lying excited state - understood as 
a systematic error - becomes negligible within the desired precision. 

Let us begin by discussing in what ways the excited states of the nu- 
cleon contribute to the nucleon matrix element calculated from the ratio 
of the 3-point to the 2-point function. The expressions given in Eqs. ([I]) 
and ([3]) illustrate how calculations of Euclidean time correlation functions 
can be used to determine matrix elements in the limit of large Euclidean 
time separations. However, in practice the finite time extension of the lat- 
tice prevents us from taking the asymptotic limits and therefore one has to 
carefully examine the sub-leading contributions usually ignored. In the fol- 
lowing we assume that at the values of t we use, excited state contributions 
in the 2-point function can be ignored. This is plausible since excited state 
contributions to the 2-point functions are generically more suppressed than 
those contributions to the 3-point function. The reasoning for this is that 
the fields in the 2-point function are always separated by a distance t and 
those in the 3-point function are separated by t — if or t' , both of which are 
smaller than t in practice. In particular, in the 3-point function we have a 
double limit requiring both t — t' and t' to be asymptotic. Let us now con- 
sider the leading contributions to the ratio of the 3-point and the 2-point 
functions originating from taking into account the contribution of the first 
excited state: 



C 3 (t,t') _ 



\O\0) 



+ (0|O|l)^y exp {-AMt') 

J AT 



>N 
7(1) 



<l|O|0> ^ exp [-AAf (t - 0] 

jv 



+ O [exp (—AMt)] , 

where AM is the mass gap between the nucleon ground state and the first 
excited state. As can be seen, there are two additional time dependent 
contributions to leading order. 
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3. Lattice techniques and details 

3.1. Lattice fermion action 

For this work, we employ maximally- twisted-mass Wilson fermions |10| . 
We use the gauge field configurations generated by the European Twisted 
Mass Collaboration (ETMC) with Nf = 2 + 1 + 1 flavors, thus fully ac- 
counting for the first two quark generations. We refer to Ref. for the 
details of our lattice formulation. Since we are aiming at a precise result, 
we concentrate on only one ensemble with a pion mass of m n ~ 380 MeV 
and a lattice spacing of a ~ 0.078 fm. This pion mass is chosen small 
enough to be relatively close to the physical pion mass but still large enough 
to ensure that finite size effects are suppressed. In addition, the propaga- 
tors can be calculated with moderate computational cost which allows us 
to analyze a large statistical ensemble in order to obtain an accurate result. 
We emphasize that maximally-twisted-mass fermions realize an automatic 
0(a)-improvement for which no additional operator specific improvements 
are needed. Therefore, at the lattice spacing employed here, one expects that 
discretization effects are also suppressed. This is confirmed by direct calcu- 
lations of these matrix elements at three different lattice spacings smaller 
than 0.1 fm in Refs. E1H2I- 

Although in this analysis we use twisted-mass fermions, the most impor- 
tant aspects of excited state contributions are expected to be universal and 
independent of the particular lattice discretization used. Thus the conclu- 
sions obtained here are of direct relevance to the calculations of other groups 
including those using different lattice actions. 

3.2. Fixed Sink Method 

An efficient computation of the connected piece of the 3-point function 
is possible by means of the so-called sequential source method. This tech- 
nique requires two sets of propagators. The first are the forward quark 
propagators of the light flavors that are also used to compute the 2-point 
function and are independent of the operator and hadronic state as long 
as the same source is used. Those propagators are then used to build a 
sequential source for a second generalized propagator, again for each of the 
light flavors, that is specific for the hadron state we are interested in. The 
desired matrix element is obtained by contracting the corresponding opera- 
tor with the free ends of these two propagators as illustrated in Fig. [TJ This 
method has the advantage that we can use the same set of propagators for 
any choice of the operator insertion and hence it is the method of choice for 
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Figure 1: Diagrammatic illustration of the sequential method through the sink (left) and 
the open sink method (right). 

the calculation of generalized form factors of a particular hadron. A disad- 
vantage of this method is that the source-sink separation time t in Eq. ^ 
must be fixed before the sequential propagator can be computed. Obviously, 
changing t would require another set of propagator computations. Thus, a 
prudent choice for the value of t is mandatory. This is because, on the one 
hand, according to Eq. ([3]), a large source-sink separation is desirable for the 
suppression of excited state contributions. On the other hand, the signal- 
to-noise ratio drops exponentially fast with the source-sink separation. A 
reasonable choice is therefore a source-sink separation at which the contribu- 
tion from excited states becomes negligible compared to the statistical error. 
However, this can only be determined a posteriori and only after repeating 
the entire calculation for several values of t. Thus in practice, having chosen 
a reasonable value for t, one looks at the time dependence of the right hand 
side of Eq. (|3j) clS cl function of the operator insertion time t'. If a plateau, as 
a function of t', is observed, then it is assumed that the excited states have 
been sufficiently suppressed and the plateau value is identified as the matrix 
element of interest. However, there still remains the possibility that the 
asymptotic plateau value has not been reached. Therefore in this paper, we 
carry out a thorough investigation of the excited state contributions using 
a more appropriate approach as described in the next section. 

3.3. Open Sink Method 

In the fixed sink method the summation over x in the 3-point function 
of Eq. ([2]) is done by the sequential inversion. One can instead perform 
the summation over y through a sequential inversion. Therefore we need 
to fix the particular operator that we are interested in and also the time 
slice t' where it is inserted. The sequential source is then constructed at 
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the operator rather than the sink and hence, the sequential propagator is 
now operator-dependent but state and t independent. For the example of 
an operator that involves only one lattice site, the sequential propagator is 
given by 

SSk(s), t; x h 0;O) = J2 S&(x* f ,t; y, t') 0%{y, t') Sf a {y, t'; x% 0) 
v 

where the Roman superscripts are color indices and the Greek subscripts are 
Dirac indices. The generalization of this expression for operators that involve 
more than one lattice site and gauge links, such as derivative operators, is 
straightforward. The sequential propagator is obtained for all source sink 
separations t, thereby allowing us to study the effects of excited states. The 
open sink or fixed current method is illustrated in Fig. [T] Clearly, this 
method is not practical when a large number of matrix elements of different 
operators must be computed since for each new operator extra inversions 
must be performed. However, for our dedicated study of excited state effects 
in gA and (x) u -d it is clearly very useful. 

3.4- Observables 

As stated above, in this work, we concentrate on two relatively simple 
but nonetheless phenomenologically very relevant quantities. The first is 
the nucleon axial charge, gA, which plays an important role in the beta de- 
cay of the neutron and appears as a low energy constant in effective chiral 
Lagrangians. It has been precisely measured and it is straightforward to 
calculate in lattice QCD using the techniques described in the previous sec- 
tions. However, the values obtained from various lattice QCD calculations 
are typically 5% to 10% lower while having themselves a statistical accuracy 
of the order of 1%, see Fig. [2] for the example of our own calculations. 

The second observable is the lowest non-trivial moment of the unpo- 
larized parton distribution function in isovector flavor combination, {x) u _ d . 
This quantity is determined phenomenologically from a global analysis of 
deep inelastic scattering data, and the discrepancy between the phenomeno- 
logical and lattice values is even larger, roughly 50% to 60%, see Fig. [2] 
For the precise definitions of the corresponding 3-point functions, we refer 
to Refs. [U |6]. Though this does not affect excited state contamination, 
it is important to note that we use a non-perturbative renormalization of 
our bare matrix elements. The corresponding renormalization factors are 
calculated in the RI'MOM scheme and are matched to the MS scheme at a 
scale of (2 GeV) 2 . For more details we refer to Refs. [HI EJ- The values 
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of the renormalization constants used in this work are Za = 0.774 for the 
renormalization of the bare gA and Z/ x \ = 0.998 for the renormalization of 

( X )u-d- 

4. Results 

In order to have a reference value, we first perform a calculation of the 
nucleon axial charge qa using the sequential fixed sink method with a fixed 
source-sink separation of t = 12a ~ 0.94 fm. Gauge invariant Gaussian 
smearing of the quark fields, including APE-smeared gauge links, is used in 
order to improve the overlap with the nucleon ground state. In Fig. [2] we 
compare the value obtained for the Nf = 2 + 1 + 1 ensemble to our previous 
Nf = 2 results at various pion masses. As can be seen, the value we find for 
qa is in good agreement with the Nf = 2 results obtained at nearby pion 
masses. 

We then perform an analysis on the same Nf = 2 + 1 + 1 ensemble using 
the open sink method. The time slice of the operator insertion was fixed to 
t 1 = 9a. This was chosen to safely suppress excited state contributions from 
the source, as can be verified from the 2-point function. The result of the 
analysis using the open sink method is shown in Fig. [3| Clearly, the value of 
gA does not show any statistically significant dependence on the source-sink 
separation t. Hence, the plot demonstrates that there is no contribution 
from excited states detectable within the statistical accuracy of 2.5%. Note 
that, although t' = 9a, a plateau for gA is already reached at t = 11a. It is 
worth mentioning that, in order to reach a comparable statistical accuracy 
as the one obtained when using the fixed sink method with t = 12a with 
500 measurements, we had to perform roughly 7500 measurements when we 
take for example t = 18a. This was achieved by choosing randomly located 
source points with typically 2 sources per configuration. 

As a second benchmark quantity, we have examined (x) >. As for the 
case of gA, we first determine the value of this observable using the fixed 
sink method increasing the statistics of the calculation presented in |14j . 
For the open sink method, we have chosen the operator insertion time to 
be t' = 11a. We expect that for this choice excited state effects from the 
source are sufficiently suppressed for this operator. We perform in total 
about 23,000 measurements for (x) u _ d using randomly distributed source 
points with 5 sources per configuration. With such statistics and at a source- 
sink separation of t = 18a, we could equal the precision of the fixed sink 
method that was done with a source-sink separation of t = 12a using 1300 
measurements. In Fig. ffl we plot (x) u _ d as a function of the source-sink 
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Figure 2: In the left panel, we show the relative deviation of ETMC lattice results for 
gA from the experimental value [15]. In the right panel, we show the relative deviation of 
ETMC lattice results for (x) d from a result obtained from a phenomenological analy- 
sis |16| . The lattice values for Nf — 2 at the various pion masses are from Refs. [HI I12j. 
The filled (magenta) diamonds show the results using the Nf = 2 + 1 + 1 ensembles. 

separation t. We also indicate the value obtained from the fixed sink method 
analysis as well as the experimental result extracted from a recent global 
analysis [16J. 

As can be seen, for this observable there is a shift of the value of {x) u _ d 
and a plateau is reached at larger values of the source-sink separation than 
what we have used in the fixed sink method. Despite the fact that the results 
for larger values of t decrease they clearly do not reach the phenomenological 
value. In order to estimate any residual dependence on t, we determined the 
value of {x) u _ d for an infinite source-sink separation by fitting the expected 
exponential behavior, 

(s>„_ d + Aexp[-AAT(t-i / )] , 

to the lattice results with a fixed t' = 11a. The result of this fit is {x) u _ d = 
0.22(1), which is 12% lower than the result of {x) u _ d = 0.250(6), obtained 
using t = 12a in the fixed sink method. The error of the fit is estimated by 
varying the fit range and by comparing the use of a fixed parameter AM 
versus fitting AM directly. 

5. Summary and Conclusions 

In this letter we have performed precision calculations of qa and (x) u _ d 
for a single ensemble of gauge field configurations with Nf = 2 + 1 + 1 
dynamical fermions employing a non-perturbative renormalization. We have 
investigated the behavior of these benchmark quantities as a function of the 
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Figure 3: Results for qa for a range of source-sink separations obtained from the open 
sink analysis on one Nf = 2 + 1 + 1 ensemble. The light grey band indicates the result 
obtained from the fixed sink method using a source-sink separation of 12a and the dark 
grey band shows the experimental value. 

source-sink separation in order to assess the influence of excited states on 
the current lattice results for qa and (x) u _ d . This is particularly important 
given that excited states may play a role in explaining the presently observed 
discrepancy between lattice computations and phenomenological evaluations 
of several important nucleon observables. 

We find that for the here considered pion mass of about 380 MeV and 
lattice spacing of a ~ 0.078 fm, the contamination of excited states is neg- 
ligible for gA, but for {x) u _ d , the effect is of the order of 10% compared 
to our previous calculations, where the source-sink separation has been set 
to about 1 fm. This is an effect larger than the finite volume and lattice 
spacing effects we observe at this value of the pion mass, volume and lat- 
tice spacing. Moreover, this demonstrates that contributions from excited 
states are operator dependent and should be investigated separately for each 
operator. 

One way to better control excited state effects is to use a variational 
method such as the generalized eigenvalue method [17} 118] . Recently, a 
new approach to deal with excited state contamination of hadronic matrix 
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Figure 4: Results for (x) u _ d for a range of source-sink separations obtained by means of 
the open sink method. The operator insertion was at a temporal separation from the source 
of t' — 11a. The value (including errors) obtained from the fixed sink method using a 
source-sink separation of 12a is indicated by the light grey band. The phenomenologically 
extracted value is shown with the dark grey band. The blue solid line corresponds to a fit 
described in the text. 

elements has been developed and applied for the B*Btt coupling in ref. |19j . 
Whether the generalized eigenvalue method can improve the calculation of 
matrix elements of the nucleon needs still to be tested, though. 

However, if the 10% shift for {x) u _ d as we found here persists at smaller 
pion masses, excited state effects can not be the single dominating system- 
atic effect responsible for the tension between lattice and phenomenology. 
Of course, we cannot exclude that at smaller values of the pion mass excited 
state effects might become significantly larger. Therefore, in order to clarify 
the deviation between lattice calculations and experimental determinations 
of nucleon matrix elements, a very careful and accurate analysis of system- 
atic errors will be needed, taking into account the possible contamination of 
excited states as observed in this work. 
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